R Notebook

# load libraries
# load libraries
library(clusterSim) 

Loading required package: cluster
Loading required package: MASS
Warning: package ‘MASS’ was built under R version 4.1.2

library(FlowSOM)

Loading required package: igraph
Warning: package ‘igraph’ was built under R version 4.1.2

Attaching package: ‘igraph’

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

Thanks for using FlowSOM. From version 2.1.4 on, the scale 
parameter in the FlowSOM function defaults to FALSE

library(flowCore)

Attaching package: ‘flowCore’

The following object is masked from ‘package:igraph’:

    normalize

library(cluster)
library(fpc)
library(clv)

Loading required package: class
Warning: package ‘class’ was built under R version 4.1.2

Attaching package: ‘class’

The following object is masked from ‘package:FlowSOM’:

    SOM

The following object is masked from ‘package:igraph’:

    knn

library(Seurat)

Warning: package ‘Seurat’ was built under R version 4.1.2
Attaching SeuratObject

library(dplyr)

Warning: package ‘dplyr’ was built under R version 4.1.2

Attaching package: ‘dplyr’

The following object is masked from ‘package:flowCore’:

    filter

The following objects are masked from ‘package:igraph’:

    as_data_frame, groups, union

The following object is masked from ‘package:MASS’:

    select

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)

Need help getting started? Try the R Graphics Cookbook: https://r-graphics.org

library(clustree)

Loading required package: ggraph

library(reshape2)

Visialize correlations on UMAP

# read in correlation

# read in Seurat object to make UMAP
seu <- readRDS("/Users/rhalenathomas/Documents/Data/FlowCytometry/PhenoID/Analysis/9MBO/prepro_outsjan20-9000cells/Figure3/cluster_parameters/LouvainFinal/FlowAlignRetroTransLouvainSeuratObject60.Rds")


# keep feature order the same as in other figures

AB <- c("CD24","CD56","CD29","CD15","CD184","CD133","CD71","CD44","GLAST","AQP4","HepaCAM", "CD140a","O4")

Add the cell labels into seurat object

unique(seu$cor.labels)

 [1] "OPC"                    "RG"                     "Astrocyte"              "Endothelial"            "Neuron-OPC"            
 [6] "OPC-Epithelial"         "Astrocyte-RG"           "StemCell"               "Epithelial"             "OPC-Neuron"            
[11] "Oligo"                  "RG-Endothelial"         "RG-Astrocyte"           "StemCell-OPC"           "Neuron"                
[16] "Epithelial-Endothelial" "StemCell-Neuron"        "Endothelial-RG"         "OPC-StemCell"           "Astrocyte-Endothelial" 
[21] "StemCell-Oligo"         "NPC"                    "Endothelial-Astrocyte"  "Epithelial-OPC"         "Endothelial-Epithelial"
[26] "Neuron-NPC"             "Epithelial-StemCell"    "RG-OPC"                 "Astrocyte-NPC"          "StemCell-Endothelial"  
[31] "OPC-Oligo"              "RG-Oligo"               "Oligo-StemCell"         "NPC-StemCell"           "Neuron-RG"             
[36] "OPC-RG"                 "NPC-Astrocyte"          "RG-NPC"                 "Epithelial-RG"          "RG-Neuron"             
[41] "RG-Epithelial"          "NPC-Neuron"             "Oligo-RG"               "Astrocyte-OPC"          "OPC-Astrocyte"         
[46] "OPC-Endothelial"        "Neuron-Endothelial"     "Neuron-Astrocyte"       "Neuron-StemCell"        "StemCell-Epithelial"   
[51] "Epithelial-NPC"         "StemCell-NPC"           "Endothelial-StemCell"   "Epithelial-Astrocyte"   "Oligo-OPC"             
[56] "Oligo-Neuron"           "Oligo-NPC"              "Epithelial-Oligo"       "OPC-NPC"                "StemCell-RG"           
[61] "unknown"                "NPC-RG"                 "Oligo-Epithelial"       "Endothelial-Neuron"     "Astrocyte-Epithelial"  
[66] "NPC-OPC"                "Astrocyte-Oligo"        "Neuron-Oligo"           "Endothelial-Oligo"      "NPC-Endothelial"       
[71] "NPC-Oligo"              "Oligo-Astrocyte"        "Endothelial-NPC"        "StemCell-Astrocyte"     "Endothelial-OPC"       
[76] "Astrocyte-StemCell"     "Oligo-Endothelial"      "Astrocyte-Neuron"       "RG-StemCell"            "NPC-Epithelial"        
[81] "Epithelial-Neuron"      "Neuron-Epithelial"     

Make UMAP

Save UMAPS

pdf(paste(outpath,"UMAPCorLables9celltypes202203018.pdf"), width = 20, height = 8)
DimPlot(seu, reduction = "umap", label = TRUE, group.by = 'cor.labels', repel = TRUE, shuffle = TRUE)

Warning: ggrepel: 25 unlabeled data points (too many overlaps). Consider increasing max.overlaps

dev.off()

null device 
          1 

Add labels

0-Epithelial 1-Astrocytes-RG 2-Neurons 3-NPC (neural stem) 4-RG-CD133 5-RG-CD184 6-Astrocytes 7-Neurons-early 8-Astrocytes 9-Glia (RG/neuron) 10-Endothelial 11-Oligodendrocytes 12-mix 13-OPC 14-Neurons-early

pdf(paste(outpath,"HeatMapcell20220318.pdf"), width = 8, height = 5)
DoHeatmap(seu, group.by = 'labels4', features = AB, size= 5,slot = "scale.data", disp.max = 2, disp.min = -1.5) + scale_fill_gradientn(colors = c("#154c79", "#eeeee4", "#e28743")) + 
  theme(axis.text.y = element_text(size = 15))

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

dev.off()

null device 
          1 

Visualize split UMAPs by batch

pdf(paste(outpath,"UMAPsplitBatch_clusterlabels2.pdf"), width = 30, height = 5)
  DimPlot(seu, split.by = 'Batch', label = TRUE)
dev.off()

null device 
          1 

Plot the correlations with less mix-cell types

Try to calculate how many of each cell type are in each cluster

# total proportion in each cell type
table(seu$RNA_snn_res.0.4)

    0     1     2     3     4     5     6     7     8     9    10    11    12    13    14 
20874 10783  8565  6789  6024  5703  3296  3033  2632  2535  2379   430   245   185   105 

table(seu$labels4)

   Epithelial  Astrocytes-1       Neurons           NPC          RG-1          RG-2  Astrocytes-2 
        20874         13415          8565          6789          6024          5703          5831 
Neurons-early   Endothelial         Oligo           mix           OPC 
         3138          2379           430           245           185 

#proportions
prop.table(table(seu$labels4))

   Epithelial  Astrocytes-1       Neurons           NPC          RG-1          RG-2  Astrocytes-2 
  0.283698932   0.182323521   0.116407078   0.092269428   0.081872299   0.077509582   0.079249232 
Neurons-early   Endothelial         Oligo           mix           OPC 
  0.042648618   0.032333034   0.005844138   0.003329800   0.002514339 

Proportions of each sample:

Use correlation labels to help label clusters

cor.lables <- as.data.frame(table(seu$RNA_snn_res.1.2, seu$cor.labels))
cor.lables$Freq <- as.double(cor.lables$Freq)
# try to make a heatmap

# not working
hm <- ggplot(data = pr.cor.lables, mapping = aes(x= Var1, y = Var2, fill = Freq))
hm

# try bar chart
ggplot(cor.lables, aes(y = Freq, x = Var1, fill = Var2)) + geom_bar(position = "stack", stat= "identity")

See the same format of UMAP

DimPlot(seu, group.by = "RNA_snn_res.1.2", label = TRUE)

Make a lot of rigde plots

RidgePlot(seu, group.by = "RNA_snn_res.1.2", features = "CD140a", log = TRUE)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.131
Warning: Removed 44076 rows containing non-finite values (stat_density_ridges).

Add some new labels from resolution 1.2

See the expression in each group in ridge plot

RidgePlot(seu, group.by = "labels4", features = AB, log = TRUE)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.105
Warning: Removed 21734 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0805
Warning: Removed 618 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0772
Warning: Removed 13691 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0901
Warning: Removed 14795 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0884
Warning: Removed 11826 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0756
Warning: Removed 15435 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0806
Warning: Removed 17617 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.114
Warning: Removed 1308 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.13
Warning: Removed 27822 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.102
Warning: Removed 21525 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.094
Warning: Removed 28052 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.102
Warning: Removed 44076 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0867
Warning: Removed 11933 rows containing non-finite values (stat_density_ridges).

RidgePlot(seu, group.by = "labels6", features = AB, log = TRUE)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.102
Warning: Removed 21734 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0764
Warning: Removed 618 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0799
Warning: Removed 13691 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0917
Warning: Removed 14795 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0787
Warning: Removed 11826 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0809
Warning: Removed 15435 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0869
Warning: Removed 17617 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.105
Warning: Removed 1308 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.119
Warning: Removed 27822 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.102
Warning: Removed 21525 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0885
Warning: Removed 28052 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.101
Warning: Removed 44076 rows containing non-finite values (stat_density_ridges).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0747
Warning: Removed 11933 rows containing non-finite values (stat_density_ridges).

Get mean expression per group/cell type of each antibody

Label transfer using Seurat

predictions <- TransferData(anchorset = anchors, refdata = seu.r$labels6, dims = 1:10)

Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels

Lets see the transfer labels

We can see how many match if we had labels already

I can see how many match between two other labelss!

table(seu.r$cor.cell.match)

FALSE  TRUE 
71242  2336 

For some reason the dotplot doesn’t work - I’ll remake the dotplot with ggplot

express.by.cluster <- as.data.frame(AverageExpression(seu.q, features = AB, group.by = 'labels', slot = 'scale.data'))
express.by.cluster <- scale(express.by.cluster)

# reformat
longData<- melt(express.by.cluster)

# heatmap

ggplot(longData, aes(x = Var2, y = Var1)) + 
  geom_raster(aes(fill=value)) + 
  scale_fill_gradient(low="grey90", high="red", na.value = "grey") +
  labs(x="Antibodies", y="Cell Types", title="Matrix") +
  theme_bw() + theme(axis.text.x=element_text(size=12, angle=90, vjust=0.3),
                     axis.text.y=element_text(size=12),
                     plot.title=element_text(size=12))

ggplot(data = df.plot, aes(x = Var1, y = Var2, color = , size = GeneRatio)) +
  geom_point() +
  scale_color_gradient(low = "red", high = "blue") 

Error in is_call(uq_expr, c("after_stat", "after_scale")) : 
  argument "uq_expr" is missing, with no default

Visualize Random Forest predictions

DimPlot(seu, group.by = 'rf.pred.march15', label = TRUE)
DimPlot(seu, group.by = 'rf.pred.march25', label = TRUE)

DimPlot(seu, group.by = 'rf.pred.march15.2', label = TRUE)

DimPlot(seu, group.by = 'rf.pred.march25.2', label = TRUE)

DimPlot(seu, group.by = 'labels6', label = TRUE)

DimPlot(seu, group.by = 'labels4', label = TRUE)